<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">

<!-- 
	Copyright (C) 2007, 2008, 2009, 2010, 2011. PARP Research Group.
	<http://perception.inf.um.es>
	University of Murcia, Spain.

	This file is part of the QVision library.

	QVision is free software: you can redistribute it and/or modify
	it under the terms of the GNU Lesser General Public License as
	published by the Free Software Foundation, version 3 of the License.

	QVision is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU Lesser General Public License for more details.

	You should have received a copy of the GNU Lesser General Public
	License along with QVision. If not, see <http://www.gnu.org/licenses/>.
-->

<html><head><meta http-equiv="content-Type" content="text/html;charset=UTF-8">
<title>QVision: Qt&#39;s Image, Video and Computer Vision Library</title>
<meta name="title" content="QVision" />
<meta name="dc.title" content="QVision" />
<meta name="url" content="http://perception.inf.um.es/QVision" />
<meta name="author" content="PARP Research Group - http://perception.inf.um.es" />
<meta name="revisit-after" content="30 DAYS"/>
<meta name="robots" content="index,follow"/>
<meta name="classification" content="*">
<meta name="rating" content="Safe For Kids">
<meta name="distribution" content="GLOBAL"/>
<meta name="description" content="Qt's Image, Video and Computer Vision Library"/>
<meta name="page-topic" content="Computer Vision research and prototype programming"/>
<meta name="geo.country" content="ES" />

<!--
Keywords:
By license:		open source, gnu, lgpl, gpl, free
By theme:		computer vision, image processing, robotics, programming, source, development
By usage:		library, toolkit, framework, prototype, application
By programming specs:	object oriented, c++, block programming, reusability, gui, graphical, parallel computing, high performance, GPU, prototyping
Interoperability with:	Qt, GSL, GNU Scientific library, OpenCV, CGAL, QWT, CUDA, mplayer, IPP, Intel Image Performance Primitives, blas, lapack
Functionallity:		image features, matrix algebra, projective geometry, mser, function minimization, function optimization, canny operator, harris operator, corner detection, performance evaluation, cpu usage, graphical interface
Main data-types:	matrix, vector, tensor, quaternion, image, polyline
Video sources:		webcam, camera, stream
Devices:		embedded, desktop computer, laptop, mini-laptop
Authors:		PARP research group. University of Murcia, Spain.
-->

<meta name="keywords" content="augmented reality, sfm, structure from motion, open source, gnu, lgpl, gpl, free, computer vision, image processing, robotics, programming, source, development, library, toolkit, framework, prototype, application, object oriented, c++, block programming, reusability, gui, graphical, parallel computing, high performance, GPU, prototyping, Qt, GSL, GNU Scientific library, OpenCV, CGAL, QWT, CUDA, mplayer, IPP, Intel Image Performance Primitives, blas, lapack, image features, matrix algebra, projective geometry, mser, function minimization, function optimization, canny operator, harris operator, corner detection, performance evaluation, cpu usage, graphical interface, matrix, vector, tensor, quaternion, image, polyline, webcam, camera, stream, embedded, desktop computer, laptop, mini-laptop, University of Murcia, Spain, PARP research group, vision por computador"/>
<meta http-equiv="keywords" content="augmented reality, sfm, structure from motion, open source, gnu, lgpl, gpl, free, computer vision, image processing, robotics, programming, source, development, library, toolkit, framework, prototype, application, object oriented, c++, block programming, reusability, gui, graphical, parallel computing, high performance, GPU, prototyping, Qt, GSL, GNU Scientific library, OpenCV, CGAL, QWT, CUDA, mplayer, IPP, Intel Image Performance Primitives, blas, lapack, image features, matrix algebra, projective geometry, mser, function minimization, function optimization, canny operator, harris operator, corner detection, performance evaluation, cpu usage, graphical interface, matrix, vector, tensor, quaternion, image, polyline, webcam, camera, stream, embedded, desktop computer, laptop, mini-laptop, University of Murcia, Spain, PARP research group, vision por computador"/>
<meta http-equiv="pragma" content="no-cache"/>
<meta http-equiv="title" content="QVision"/>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="tabs.css" rel="stylesheet" type="text/css" />
<link rel="shortcut icon" href="favicon.ico" />
</head><body>

<table width="100%"><tr>
	<td><a href="http://perception.inf.um.es/"><img src="parp.png" border="0" /> <big>PARP Research Group</big></a></td>
	<td align="right"><a href="http://www.um.es/"><big>Universidad de Murcia</big> <img src="um.png" border="0" /></a></td>
</tr></table>

<hr /><br />

<table width="95%" align="center"><tr><td>

<!-- Generated by Doxygen 1.6.3 -->
<script type="text/javascript"><!--
var searchBox = new SearchBox("searchBox", "search",false,'Search');
--></script>
<div class="navigation" id="top">
  <div class="tabs">
    <ul>
      <li><a href="index.html"><span>Main&nbsp;Page</span></a></li>
      <li><a href="pages.html"><span>Related&nbsp;Pages</span></a></li>
      <li><a href="modules.html"><span>Modules</span></a></li>
      <li><a href="annotated.html"><span>Classes</span></a></li>
      <li class="current"><a href="files.html"><span>Files</span></a></li>
      <li>
        <div id="MSearchBox" class="MSearchBoxInactive">
          <form id="FSearchBox" action="search.php" method="get">
            <img id="MSearchSelect" src="search/search.png" alt=""/>
            <input type="text" id="MSearchField" name="query" value="Search" size="20" accesskey="S" 
                   onfocus="searchBox.OnSearchFieldFocus(true)" 
                   onblur="searchBox.OnSearchFieldFocus(false)"/>
          </form>
        </div>
      </li>
    </ul>
  </div>
  <div class="tabs">
    <ul>
      <li><a href="files.html"><span>File&nbsp;List</span></a></li>
      <li><a href="globals.html"><span>File&nbsp;Members</span></a></li>
    </ul>
  </div>
<h1>src/qvmath/qvnumericalanalysis.cpp</h1><a href="qvnumericalanalysis_8cpp.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/*</span>
<a name="l00002"></a>00002 <span class="comment"> *      Copyright (C) 2007, 2008, 2009, 2010, 2011, 2012. PARP Research Group.</span>
<a name="l00003"></a>00003 <span class="comment"> *      &lt;http://perception.inf.um.es&gt;</span>
<a name="l00004"></a>00004 <span class="comment"> *      University of Murcia, Spain.</span>
<a name="l00005"></a>00005 <span class="comment"> *</span>
<a name="l00006"></a>00006 <span class="comment"> *      This file is part of the QVision library.</span>
<a name="l00007"></a>00007 <span class="comment"> *</span>
<a name="l00008"></a>00008 <span class="comment"> *      QVision is free software: you can redistribute it and/or modify</span>
<a name="l00009"></a>00009 <span class="comment"> *      it under the terms of the GNU Lesser General Public License as</span>
<a name="l00010"></a>00010 <span class="comment"> *      published by the Free Software Foundation, version 3 of the License.</span>
<a name="l00011"></a>00011 <span class="comment"> *</span>
<a name="l00012"></a>00012 <span class="comment"> *      QVision is distributed in the hope that it will be useful,</span>
<a name="l00013"></a>00013 <span class="comment"> *      but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
<a name="l00014"></a>00014 <span class="comment"> *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the</span>
<a name="l00015"></a>00015 <span class="comment"> *      GNU Lesser General Public License for more details.</span>
<a name="l00016"></a>00016 <span class="comment"> *</span>
<a name="l00017"></a>00017 <span class="comment"> *      You should have received a copy of the GNU Lesser General Public</span>
<a name="l00018"></a>00018 <span class="comment"> *      License along with QVision. If not, see &lt;http://www.gnu.org/licenses/&gt;.</span>
<a name="l00019"></a>00019 <span class="comment"> */</span>
<a name="l00020"></a>00020 
<a name="l00021"></a>00021 
<a name="l00025"></a>00025 
<a name="l00026"></a>00026 <span class="preprocessor">#include &lt;<a class="code" href="qvmath_2qvnumericalanalysis_8h.html" title="File from the QVision library.">qvmath/qvnumericalanalysis.h</a>&gt;</span>
<a name="l00027"></a>00027 
<a name="l00028"></a><a class="code" href="group__qvnumericalanalysis.html#ga4d2e455ca93f788e1ec97d942171c8ef">00028</a> <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> <a class="code" href="group__qvnumericalanalysis.html#ga4d2e455ca93f788e1ec97d942171c8ef" title="Estimates the gradient vector for the function using the forward two-points rule...">qvEstimateGradient</a>(<a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> &amp;multivariateFunction, <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;location, <span class="keyword">const</span> <span class="keywordtype">double</span> h)
<a name="l00029"></a>00029     {
<a name="l00030"></a>00030     <span class="keyword">const</span> <span class="keywordtype">int</span> dim = location.size();
<a name="l00031"></a>00031     <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> gradient(dim);
<a name="l00032"></a>00032     <span class="keyword">const</span> <span class="keywordtype">double</span> actual = multivariateFunction(location);
<a name="l00033"></a>00033     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt; dim; i++)
<a name="l00034"></a>00034         {
<a name="l00035"></a>00035         <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> stepLocation = location;
<a name="l00036"></a>00036         stepLocation[i] += h;
<a name="l00037"></a>00037         gradient[i] = (multivariateFunction(stepLocation) - actual)/h;
<a name="l00038"></a>00038         }
<a name="l00039"></a>00039     <span class="keywordflow">return</span> gradient;
<a name="l00040"></a>00040     }
<a name="l00041"></a>00041 
<a name="l00042"></a><a class="code" href="group__qvnumericalanalysis.html#gadd02b7275fd52c76e9cb406031b4392c">00042</a> <span class="keyword">const</span> <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a> <a class="code" href="group__qvnumericalanalysis.html#gadd02b7275fd52c76e9cb406031b4392c" title="Estimates the Jacobian matrix for the function using the forward two-points rule...">qvEstimateJacobian</a>(<a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> &amp;multivariateFunction, <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;location, <span class="keyword">const</span> <span class="keywordtype">double</span> h)
<a name="l00043"></a>00043     {
<a name="l00044"></a>00044     <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> actual = multivariateFunction(location);
<a name="l00045"></a>00045 
<a name="l00046"></a>00046     <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a> jacobian(actual.size(), location.size());
<a name="l00047"></a>00047 
<a name="l00048"></a>00048     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt; location.size(); i++)
<a name="l00049"></a>00049         {
<a name="l00050"></a>00050         <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> stepLocation = location;
<a name="l00051"></a>00051         stepLocation[i] += h;
<a name="l00052"></a>00052         jacobian.setCol(i, (multivariateFunction(stepLocation) - actual)/h);
<a name="l00053"></a>00053         }
<a name="l00054"></a>00054 
<a name="l00055"></a>00055     <span class="keywordflow">return</span> jacobian;
<a name="l00056"></a>00056     }
<a name="l00057"></a>00057 
<a name="l00058"></a><a class="code" href="group__qvnumericalanalysis.html#gabe23aebe828a863a2d5bc305b5258179">00058</a> <span class="keyword">const</span> <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a> <a class="code" href="group__qvnumericalanalysis.html#gabe23aebe828a863a2d5bc305b5258179" title="Estimates the hessian matrix for the function using the forward two-point rule for...">qvEstimateHessian</a>(       <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> &amp;multivariateFunction,
<a name="l00059"></a>00059                     <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;location, <span class="keyword">const</span> <span class="keywordtype">double</span> h)
<a name="l00060"></a>00060     {
<a name="l00061"></a>00061     <span class="keyword">const</span> <span class="keywordtype">int</span> dim = location.size();
<a name="l00062"></a>00062     <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> g = <a class="code" href="group__qvnumericalanalysis.html#ga4d2e455ca93f788e1ec97d942171c8ef" title="Estimates the gradient vector for the function using the forward two-points rule...">qvEstimateGradient</a>(multivariateFunction, location, h);
<a name="l00063"></a>00063 
<a name="l00064"></a>00064     <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a> hessian(dim, dim);
<a name="l00065"></a>00065 
<a name="l00066"></a>00066     <span class="keyword">const</span> <span class="keywordtype">double</span> actual = multivariateFunction(location);
<a name="l00067"></a>00067     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt; dim; i++)
<a name="l00068"></a>00068         <span class="keywordflow">for</span> (<span class="keywordtype">int</span> j = 0; j &lt; dim; j++)
<a name="l00069"></a>00069         {
<a name="l00070"></a>00070         <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> stepLocationIJ = location;
<a name="l00071"></a>00071         stepLocationIJ[i] += h;
<a name="l00072"></a>00072         stepLocationIJ[j] += h;
<a name="l00073"></a>00073         hessian(i,j) = (multivariateFunction(stepLocationIJ) - actual)/(h*h) - (g[i] + g[j])/h;
<a name="l00074"></a>00074         }
<a name="l00075"></a>00075     <span class="keywordflow">return</span> hessian;
<a name="l00076"></a>00076     }
<a name="l00077"></a>00077 
<a name="l00078"></a>00078 <span class="preprocessor">#ifdef GSL_AVAILABLE</span>
<a name="l00079"></a>00079 <span class="preprocessor"></span>
<a name="l00080"></a>00080 <span class="comment">// GSL minimization</span>
<a name="l00081"></a>00081 
<a name="l00082"></a>00082 <span class="keywordtype">double</span> my_f (<span class="keyword">const</span> gsl_vector *v, <span class="keywordtype">void</span> *params)
<a name="l00083"></a>00083     {
<a name="l00084"></a>00084     <span class="keywordflow">return</span> ((<a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> *) params)-&gt;operator()(<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(v));
<a name="l00085"></a>00085     }
<a name="l00086"></a>00086 
<a name="l00087"></a>00087 <span class="comment">/* The gradient of f, df = (df/dx, df/dy). */</span>
<a name="l00088"></a>00088 <span class="keywordtype">void</span> my_df (<span class="keyword">const</span> gsl_vector *v, <span class="keywordtype">void</span> *params, gsl_vector *df)
<a name="l00089"></a>00089     {
<a name="l00090"></a>00090     gsl_vector_memcpy(df, <a class="code" href="group__qvnumericalanalysis.html#ga4d2e455ca93f788e1ec97d942171c8ef" title="Estimates the gradient vector for the function using the forward two-points rule...">qvEstimateGradient</a>( * (<a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> *) params,<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(v)));
<a name="l00091"></a>00091     }
<a name="l00092"></a>00092 
<a name="l00093"></a>00093 <span class="comment">/* Compute both f and df together. */</span>
<a name="l00094"></a>00094 <span class="keywordtype">void</span> my_fdf (<span class="keyword">const</span> gsl_vector *x, <span class="keywordtype">void</span> *params, <span class="keywordtype">double</span> *f, gsl_vector *df)
<a name="l00095"></a>00095     {
<a name="l00096"></a>00096     *f = my_f(x, params);
<a name="l00097"></a>00097     my_df(x, params, df);
<a name="l00098"></a>00098     }
<a name="l00099"></a>00099 
<a name="l00100"></a><a class="code" href="group__qvnumericalanalysis.html#gad5edc708c8725e70dc43ad236729d6e5">00100</a> <span class="keywordtype">bool</span> <a class="code" href="group__qvnumericalanalysis.html#gad5edc708c8725e70dc43ad236729d6e5" title="Wrapper to GSL multivariate function minimization using gradient information.">qvGSLMinimizeFDF</a> ( <span class="keyword">const</span> <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> &amp; function, <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;point,
<a name="l00101"></a>00101                 <span class="keyword">const</span> <a class="code" href="group__qvnumericalanalysis.html#gad24495d3074466956a7da56fabf2d0e7" title="GSL multidimensional minimization algorithms using gradient information.">GSLMultiminFDFMinimizerType</a> gslMinimizerAlgorithm,
<a name="l00102"></a>00102                 <span class="keyword">const</span> <span class="keywordtype">int</span> maxIterations, <span class="keyword">const</span> <span class="keywordtype">double</span> maxGradientNorm,
<a name="l00103"></a>00103                 <span class="keyword">const</span> <span class="keywordtype">double</span> step, <span class="keyword">const</span> <span class="keywordtype">double</span> tol)
<a name="l00104"></a>00104     {
<a name="l00105"></a>00105     <span class="keyword">const</span> <span class="keywordtype">int</span> dims = point.size();
<a name="l00106"></a>00106     <span class="keyword">const</span> gsl_multimin_fdfminimizer_type *minimizer_type = NULL;
<a name="l00107"></a>00107     <span class="keywordflow">switch</span>(gslMinimizerAlgorithm)
<a name="l00108"></a>00108         {
<a name="l00109"></a>00109         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a5d65b522cf13b9975a0fee9ab4f07b0e" title="Fletcher-Reeves conjugate gradient algorithm.">ConjugateFR</a>:       minimizer_type = gsl_multimin_fdfminimizer_conjugate_fr;        <span class="keywordflow">break</span>;
<a name="l00110"></a>00110         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a2c6a593f8b4eafd33170df4f5db2bc9d" title="Polak-Ribiere conjugate gradient algorithm.">ConjugatePR</a>:       minimizer_type = gsl_multimin_fdfminimizer_conjugate_pr;        <span class="keywordflow">break</span>;
<a name="l00111"></a>00111         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a1af5c2e022b228fc7e5b331e6d335c4b" title="Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.">VectorBFGS</a>:        minimizer_type = gsl_multimin_fdfminimizer_vector_bfgs;         <span class="keywordflow">break</span>;
<a name="l00112"></a>00112         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a815b31aa2c26e45abdbaf05ef58c1943" title="The steepest descent algorithm.">SteepestDescent</a>:   minimizer_type = gsl_multimin_fdfminimizer_steepest_descent;    <span class="keywordflow">break</span>;
<a name="l00113"></a>00113         }
<a name="l00114"></a>00114 
<a name="l00115"></a>00115     gsl_multimin_fdfminimizer *minimizer = gsl_multimin_fdfminimizer_alloc (minimizer_type, dims);
<a name="l00116"></a>00116 
<a name="l00117"></a>00117     gsl_multimin_function_fdf my_func;
<a name="l00118"></a>00118     my_func.n = dims;
<a name="l00119"></a>00119     my_func.f = &amp;my_f;
<a name="l00120"></a>00120     my_func.df = &amp;my_df;
<a name="l00121"></a>00121     my_func.fdf = &amp;my_fdf;
<a name="l00122"></a>00122     my_func.params = <span class="keyword">const_cast&lt;</span><a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> *<span class="keyword">&gt;</span>(&amp;function);
<a name="l00123"></a>00123 
<a name="l00124"></a>00124     gsl_vector *x = point;
<a name="l00125"></a>00125 
<a name="l00126"></a>00126     gsl_multimin_fdfminimizer_set (minimizer, &amp;my_func, x, step, tol);
<a name="l00127"></a>00127 
<a name="l00128"></a>00128     <span class="keywordtype">int</span> status = GSL_CONTINUE;
<a name="l00129"></a>00129     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; status == GSL_CONTINUE &amp;&amp; i &lt; maxIterations; i++)
<a name="l00130"></a>00130         {
<a name="l00131"></a>00131         <span class="keywordflow">if</span> ((status = gsl_multimin_fdfminimizer_iterate (minimizer)))
<a name="l00132"></a>00132             <span class="keywordflow">break</span>;
<a name="l00133"></a>00133 
<a name="l00134"></a>00134         status = gsl_multimin_test_gradient (gsl_multimin_fdfminimizer_gradient(minimizer), maxGradientNorm);
<a name="l00135"></a>00135         }
<a name="l00136"></a>00136 
<a name="l00137"></a>00137     <span class="comment">// Store resulting value.</span>
<a name="l00138"></a>00138     point = <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(gsl_multimin_fdfminimizer_x(minimizer));
<a name="l00139"></a>00139 
<a name="l00140"></a>00140     gsl_multimin_fdfminimizer_free (minimizer);
<a name="l00141"></a>00141 
<a name="l00142"></a>00142     gsl_vector_free (x);
<a name="l00143"></a>00143 
<a name="l00144"></a>00144     <span class="keywordflow">return</span> (status == GSL_SUCCESS);
<a name="l00145"></a>00145     }
<a name="l00146"></a>00146 
<a name="l00148"></a>00148 
<a name="l00149"></a>00149 <span class="keywordtype">double</span> my_f_gradient (<span class="keyword">const</span> gsl_vector *v, <span class="keywordtype">void</span> *params)
<a name="l00150"></a>00150     {
<a name="l00151"></a>00151     <span class="keywordflow">return</span> ((QPair&lt; <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> *, <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> *&gt; *) params)-&gt;first-&gt;operator()(<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(v));
<a name="l00152"></a>00152     }
<a name="l00153"></a>00153 
<a name="l00154"></a>00154 <span class="comment">/* The gradient of f, df = (df/dx, df/dy). */</span>
<a name="l00155"></a>00155 <span class="keywordtype">void</span> my_df_gradient (<span class="keyword">const</span> gsl_vector *v, <span class="keywordtype">void</span> *params, gsl_vector *df)
<a name="l00156"></a>00156     {
<a name="l00157"></a>00157     gsl_vector_memcpy(df, ((QPair&lt; <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> *, <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> *&gt; *) params)-&gt;second-&gt;operator()(<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(v)));
<a name="l00158"></a>00158     }
<a name="l00159"></a>00159 
<a name="l00160"></a>00160 <span class="comment">/* Compute both f and df together. */</span>
<a name="l00161"></a>00161 <span class="keywordtype">void</span> my_fdf_gradient (<span class="keyword">const</span> gsl_vector *x, <span class="keywordtype">void</span> *params, <span class="keywordtype">double</span> *f, gsl_vector *df)
<a name="l00162"></a>00162     {
<a name="l00163"></a>00163     *f = my_f_gradient(x, params);
<a name="l00164"></a>00164     my_df_gradient(x, params, df);
<a name="l00165"></a>00165     }
<a name="l00166"></a>00166 
<a name="l00167"></a><a class="code" href="group__qvnumericalanalysis.html#ga0dfcc84308d4ec8b2cb2d11a59251e9c">00167</a> <span class="keywordtype">bool</span> <a class="code" href="group__qvnumericalanalysis.html#gad5edc708c8725e70dc43ad236729d6e5" title="Wrapper to GSL multivariate function minimization using gradient information.">qvGSLMinimizeFDF</a> ( <span class="keyword">const</span> <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;QVVector, double&gt;</a> &amp; function, <span class="keyword">const</span> <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> &amp; gradientFunction,
<a name="l00168"></a>00168                 <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;point, <span class="keyword">const</span> <a class="code" href="group__qvnumericalanalysis.html#gad24495d3074466956a7da56fabf2d0e7" title="GSL multidimensional minimization algorithms using gradient information.">GSLMultiminFDFMinimizerType</a> gslMinimizerAlgorithm,
<a name="l00169"></a>00169                 <span class="keyword">const</span> <span class="keywordtype">int</span> maxIterations, <span class="keyword">const</span> <span class="keywordtype">double</span> maxGradientNorm,
<a name="l00170"></a>00170                 <span class="keyword">const</span> <span class="keywordtype">double</span> step, <span class="keyword">const</span> <span class="keywordtype">double</span> tol)
<a name="l00171"></a>00171     {
<a name="l00172"></a>00172     <span class="keyword">const</span> <span class="keywordtype">int</span> dims = point.size();
<a name="l00173"></a>00173     <span class="keyword">const</span> gsl_multimin_fdfminimizer_type *minimizer_type = NULL;
<a name="l00174"></a>00174     <span class="keywordflow">switch</span>(gslMinimizerAlgorithm)
<a name="l00175"></a>00175         {
<a name="l00176"></a>00176         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a5d65b522cf13b9975a0fee9ab4f07b0e" title="Fletcher-Reeves conjugate gradient algorithm.">ConjugateFR</a>:       minimizer_type = gsl_multimin_fdfminimizer_conjugate_fr;        <span class="keywordflow">break</span>;
<a name="l00177"></a>00177         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a2c6a593f8b4eafd33170df4f5db2bc9d" title="Polak-Ribiere conjugate gradient algorithm.">ConjugatePR</a>:       minimizer_type = gsl_multimin_fdfminimizer_conjugate_pr;        <span class="keywordflow">break</span>;
<a name="l00178"></a>00178         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a1af5c2e022b228fc7e5b331e6d335c4b" title="Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.">VectorBFGS</a>:        minimizer_type = gsl_multimin_fdfminimizer_vector_bfgs;         <span class="keywordflow">break</span>;
<a name="l00179"></a>00179         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggad24495d3074466956a7da56fabf2d0e7a815b31aa2c26e45abdbaf05ef58c1943" title="The steepest descent algorithm.">SteepestDescent</a>:   minimizer_type = gsl_multimin_fdfminimizer_steepest_descent;    <span class="keywordflow">break</span>;
<a name="l00180"></a>00180         }
<a name="l00181"></a>00181 
<a name="l00182"></a>00182     gsl_multimin_fdfminimizer *minimizer = gsl_multimin_fdfminimizer_alloc (minimizer_type, dims);
<a name="l00183"></a>00183     QPair&lt; const QVFunction&lt;QVVector, double&gt; *, <span class="keyword">const</span> <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> *&gt; functions(&amp;function, &amp;gradientFunction);
<a name="l00184"></a>00184 
<a name="l00185"></a>00185     gsl_multimin_function_fdf my_func;
<a name="l00186"></a>00186     my_func.n = dims;
<a name="l00187"></a>00187     my_func.f = &amp;my_f_gradient;
<a name="l00188"></a>00188     my_func.df = &amp;my_df_gradient;
<a name="l00189"></a>00189     my_func.fdf = &amp;my_fdf_gradient;
<a name="l00190"></a>00190     my_func.params = &amp;functions;
<a name="l00191"></a>00191 
<a name="l00192"></a>00192     gsl_vector *x = point;
<a name="l00193"></a>00193 
<a name="l00194"></a>00194     gsl_multimin_fdfminimizer_set (minimizer, &amp;my_func, x, step, tol);
<a name="l00195"></a>00195 
<a name="l00196"></a>00196     <span class="keywordtype">int</span> status = GSL_CONTINUE;
<a name="l00197"></a>00197     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; status == GSL_CONTINUE &amp;&amp; i &lt; maxIterations; i++)
<a name="l00198"></a>00198         {
<a name="l00199"></a>00199         <span class="keywordflow">if</span> ((status = gsl_multimin_fdfminimizer_iterate (minimizer)))
<a name="l00200"></a>00200             <span class="keywordflow">break</span>;
<a name="l00201"></a>00201 
<a name="l00202"></a>00202         status = gsl_multimin_test_gradient (gsl_multimin_fdfminimizer_gradient(minimizer), maxGradientNorm);
<a name="l00203"></a>00203         }
<a name="l00204"></a>00204 
<a name="l00205"></a>00205     <span class="comment">// Store resulting value.</span>
<a name="l00206"></a>00206     point = <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(gsl_multimin_fdfminimizer_x(minimizer));
<a name="l00207"></a>00207 
<a name="l00208"></a>00208     gsl_multimin_fdfminimizer_free (minimizer);
<a name="l00209"></a>00209 
<a name="l00210"></a>00210     gsl_vector_free (x);
<a name="l00211"></a>00211 
<a name="l00212"></a>00212     <span class="keywordflow">return</span> (status == GSL_SUCCESS);
<a name="l00213"></a>00213     }
<a name="l00214"></a>00214 
<a name="l00215"></a>00215 <span class="keywordtype">int</span> multifit_f (<span class="keyword">const</span> gsl_vector * x, <span class="keywordtype">void</span> *params, gsl_vector * f)
<a name="l00216"></a>00216     {
<a name="l00217"></a>00217     <span class="keyword">const</span> <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> result = ((QPair&lt; QVFunction&lt;QVVector, QVVector&gt; *, <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVMatrix&gt;</a> *&gt; *)params)-&gt;first-&gt;operator()(<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(x));
<a name="l00218"></a>00218 
<a name="l00219"></a>00219     gsl_vector_memcpy(f, result);
<a name="l00220"></a>00220 
<a name="l00221"></a>00221     <span class="keywordflow">return</span> GSL_SUCCESS;
<a name="l00222"></a>00222     }
<a name="l00223"></a>00223 
<a name="l00224"></a>00224 <span class="keywordtype">int</span> multifit_df (<span class="keyword">const</span> gsl_vector * x, <span class="keywordtype">void</span> *params, gsl_matrix * J)
<a name="l00225"></a>00225     {
<a name="l00226"></a>00226     <span class="keyword">const</span> <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a> result = ((QPair&lt; QVFunction&lt;QVVector, QVVector&gt; *, <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVMatrix&gt;</a> *&gt; *)params)-&gt;second-&gt;operator()(<a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a>(x));
<a name="l00227"></a>00227 
<a name="l00228"></a>00228     gsl_matrix_memcpy(J, result);
<a name="l00229"></a>00229 
<a name="l00230"></a>00230     <span class="keywordflow">return</span> GSL_SUCCESS;
<a name="l00231"></a>00231     }
<a name="l00232"></a>00232 
<a name="l00233"></a>00233 <span class="keywordtype">int</span> multifit_fdf (<span class="keyword">const</span> gsl_vector * x, <span class="keywordtype">void</span> *params, gsl_vector * f, gsl_matrix * J)
<a name="l00234"></a>00234     {
<a name="l00235"></a>00235     multifit_f (x, params, f);
<a name="l00236"></a>00236     multifit_df (x, params, J);
<a name="l00237"></a>00237 
<a name="l00238"></a>00238     <span class="keywordflow">return</span> GSL_SUCCESS;
<a name="l00239"></a>00239     }
<a name="l00240"></a>00240 
<a name="l00241"></a>00241 <span class="keywordtype">bool</span> <a class="code" href="group__qvnumericalanalysis.html#ga035c2ca4bd03efc36129a9c618e18eee" title="Solves a non-linear system of equations.">qvGSLSolveFDF</a> (    <span class="keyword">const</span> <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVVector&gt;</a> &amp; function, <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVMatrix&gt;</a> &amp; functionJacobian,
<a name="l00242"></a>00242                 <a class="code" href="classQVVector.html" title="Implementation of numerical vectors.">QVVector</a> &amp;x, <span class="keyword">const</span> <a class="code" href="group__qvnumericalanalysis.html#gaa1d694d467425bf3cfbdf949ef0faed5" title="GSL multidimensional solving.">GSLMultiminFDFSolverType</a> gslSolverAlgorithm,
<a name="l00243"></a>00243                 <span class="keyword">const</span> <span class="keywordtype">int</span> maxIterations, <span class="keyword">const</span> <span class="keywordtype">double</span> maxAbsErr, <span class="keyword">const</span> <span class="keywordtype">double</span> maxRelErr)
<a name="l00244"></a>00244     {
<a name="l00245"></a>00245     <span class="keyword">const</span> gsl_multifit_fdfsolver_type *solver_type = NULL;
<a name="l00246"></a>00246     <span class="keywordflow">switch</span>(gslSolverAlgorithm)
<a name="l00247"></a>00247         {
<a name="l00248"></a>00248         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggaa1d694d467425bf3cfbdf949ef0faed5a1bbb22ad7e6a45514c3e6022a2e50526" title="Scaled Levenberg-Marquardt algorithm.">LMScaledDerivative</a>:        solver_type = gsl_multifit_fdfsolver_lmsder;    <span class="keywordflow">break</span>;
<a name="l00249"></a>00249         <span class="keywordflow">case</span> <a class="code" href="group__qvnumericalanalysis.html#ggaa1d694d467425bf3cfbdf949ef0faed5a47a6ed26f8d738391705fb14057484eb" title="Non-scaled (faster) Levenberg-Marquardt algorithm.">LMDerivative</a>:              solver_type = gsl_multifit_fdfsolver_lmder;     <span class="keywordflow">break</span>;
<a name="l00250"></a>00250         }
<a name="l00251"></a>00251 
<a name="l00252"></a>00252     gsl_vector *x_gsl = x;
<a name="l00253"></a>00253 
<a name="l00254"></a>00254     QPair&lt; const QVFunction&lt;QVVector, QVVector&gt; *, <span class="keyword">const</span> <a class="code" href="classQVFunction.html">QVFunction&lt;QVVector, QVMatrix&gt;</a> *&gt; functions(&amp;function, &amp;functionJacobian);
<a name="l00255"></a>00255 
<a name="l00256"></a>00256     <span class="comment">// Just to get the input and output vectors sizes, for the function.</span>
<a name="l00257"></a>00257     <span class="keyword">const</span> <a class="code" href="classQVMatrix.html" title="Implementation of numerical matrices.">QVMatrix</a>      jacobian = functionJacobian(x);
<a name="l00258"></a>00258     <span class="keyword">const</span> <span class="keywordtype">int</span>   inputVectorSize = jacobian.<a class="code" href="classQVMatrix.html#a420bba03aeccbd18161418049a025f66" title="Get width of the matrix.">getCols</a>(),   <span class="comment">// Size of input vector for function.</span>
<a name="l00259"></a>00259             outputVectorSize = jacobian.<a class="code" href="classQVMatrix.html#a4108aa685baecab8a9822dcc04e98b7f" title="Get height of the matrix.">getRows</a>();      <span class="comment">// Size of output vector for function.</span>
<a name="l00260"></a>00260 
<a name="l00261"></a>00261     <span class="comment">// Generate minimization problem.</span>
<a name="l00262"></a>00262     gsl_multifit_function_fdf f;
<a name="l00263"></a>00263     f.f = &amp;multifit_f;
<a name="l00264"></a>00264     f.df = &amp;multifit_df;
<a name="l00265"></a>00265     f.fdf = &amp;multifit_fdf;
<a name="l00266"></a>00266     f.n = outputVectorSize;
<a name="l00267"></a>00267     f.p = inputVectorSize;
<a name="l00268"></a>00268     f.params = &amp;functions;
<a name="l00269"></a>00269 
<a name="l00270"></a>00270     gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc (solver_type, outputVectorSize, inputVectorSize);
<a name="l00271"></a>00271     gsl_multifit_fdfsolver_set (s, &amp;f, x_gsl);
<a name="l00272"></a>00272 
<a name="l00273"></a>00273     <span class="comment">// Perform minimization.</span>
<a name="l00274"></a>00274     <span class="keywordtype">int</span> iter = 0, status;
<a name="l00275"></a>00275     <span class="keywordflow">do</span>  {
<a name="l00276"></a>00276         iter++;
<a name="l00277"></a>00277         status = gsl_multifit_fdfsolver_iterate (s);
<a name="l00278"></a>00278 
<a name="l00279"></a>00279         <span class="keywordflow">if</span> (status)
<a name="l00280"></a>00280             <span class="keywordflow">break</span>;
<a name="l00281"></a>00281 
<a name="l00282"></a>00282         status = gsl_multifit_test_delta (s-&gt;dx, s-&gt;x, maxAbsErr, maxRelErr);
<a name="l00283"></a>00283         }
<a name="l00284"></a>00284     <span class="keywordflow">while</span> (status == GSL_CONTINUE &amp;&amp; iter &lt; maxIterations);
<a name="l00285"></a>00285 
<a name="l00286"></a>00286     <span class="comment">// Print output data.</span>
<a name="l00287"></a>00287     x = s-&gt;x;
<a name="l00288"></a>00288 
<a name="l00289"></a>00289     gsl_multifit_fdfsolver_free (s);
<a name="l00290"></a>00290     gsl_vector_free (x_gsl);
<a name="l00291"></a>00291 
<a name="l00292"></a>00292     <span class="keywordflow">return</span> (status == GSL_SUCCESS);
<a name="l00293"></a>00293     }
<a name="l00294"></a>00294 
<a name="l00296"></a>00296 
<a name="l00297"></a>00297 <span class="keywordtype">double</span> fn2 (<span class="keywordtype">double</span> x, <span class="keywordtype">void</span> * params)
<a name="l00298"></a>00298     {
<a name="l00299"></a>00299     <span class="keywordflow">return</span> ((<a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;double, double&gt;</a> *) params)-&gt;operator()(x);
<a name="l00300"></a>00300     }
<a name="l00301"></a>00301 
<a name="l00302"></a><a class="code" href="group__qvnumericalanalysis.html#ga19129c12d25402096905b2cd34529d26">00302</a> <span class="keywordtype">bool</span> <a class="code" href="group__qvnumericalanalysis.html#ga19129c12d25402096905b2cd34529d26" title="Wrapper to GSL function minimization.">qvGSLMinimize</a>(<span class="keyword">const</span> <a class="code" href="classQVFunction.html" title="Base class for function objects.">QVFunction&lt;double, double&gt;</a> &amp;function,
<a name="l00303"></a>00303                          <span class="keywordtype">double</span> &amp;x, <span class="keywordtype">double</span> &amp;lower, <span class="keywordtype">double</span> &amp;upper,
<a name="l00304"></a>00304                          <span class="keyword">const</span> <a class="code" href="group__qvnumericalanalysis.html#gaf94804b1e1dc7c0e9af0999d68440f70" title="GSL Minimization algorithms.">GSLMinFMinimizer</a> gslMinimizerAlgorithm,
<a name="l00305"></a>00305                          <span class="keyword">const</span> <span class="keywordtype">int</span> maxIterations,
<a name="l00306"></a>00306                          <span class="keyword">const</span> <span class="keywordtype">double</span> absoluteError,
<a name="l00307"></a>00307                          <span class="keyword">const</span> <span class="keywordtype">double</span> relativeError)
<a name="l00308"></a>00308     {
<a name="l00309"></a>00309     <span class="keyword">const</span> gsl_min_fminimizer_type *minimizer_type =
<a name="l00310"></a>00310             (gslMinimizerAlgorithm == <a class="code" href="group__qvnumericalanalysis.html#ggaf94804b1e1dc7c0e9af0999d68440f70ae769900a9fc183f3501fc0fdc412ddb9" title="The golden section algorithm. The simplest method of bracketing the minimum of a...">GoldenSection</a>)? gsl_min_fminimizer_goldensection : gsl_min_fminimizer_brent;
<a name="l00311"></a>00311 
<a name="l00312"></a>00312     gsl_function F;
<a name="l00313"></a>00313     F.function = &amp;fn2;
<a name="l00314"></a>00314     F.params = (<span class="keywordtype">void</span> *) &amp;function;
<a name="l00315"></a>00315 
<a name="l00316"></a>00316     gsl_min_fminimizer *s = gsl_min_fminimizer_alloc (minimizer_type);
<a name="l00317"></a>00317     gsl_min_fminimizer_set (s, &amp;F, x, lower, upper);
<a name="l00318"></a>00318 
<a name="l00319"></a>00319     <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i = 0; i &lt;maxIterations; i++)
<a name="l00320"></a>00320         {
<a name="l00321"></a>00321         gsl_min_fminimizer_iterate(s);
<a name="l00322"></a>00322         x = gsl_min_fminimizer_x_minimum (s);
<a name="l00323"></a>00323         lower = gsl_min_fminimizer_x_lower (s);
<a name="l00324"></a>00324         upper = gsl_min_fminimizer_x_upper (s);
<a name="l00325"></a>00325         <span class="keywordflow">if</span> (gsl_min_test_interval(lower, upper, absoluteError, relativeError) != GSL_CONTINUE)
<a name="l00326"></a>00326             <span class="keywordflow">break</span>;
<a name="l00327"></a>00327         }
<a name="l00328"></a>00328 
<a name="l00329"></a>00329     gsl_min_fminimizer_free (s);
<a name="l00330"></a>00330 
<a name="l00331"></a>00331     <span class="keywordflow">return</span> gsl_min_test_interval(lower, upper, absoluteError, relativeError) == GSL_SUCCESS;
<a name="l00332"></a>00332     }
<a name="l00333"></a>00333 <span class="preprocessor">#endif</span>
</pre></div></div>
</td></tr></table>

<br /><hr><br />
<center><a href="http://perception.inf.um.es/QVision">QVision framework</a>.
<a href="http://perception.inf.um.es">PARP research group</a>.
Copyright &copy; 2007, 2008, 2009, 2010, 2011.</center>
<br />
</body>
</html>
